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ADAPTATION OF A GENERAL CIRCULATION 
MODEL TO OCEAN DYNAMICS 

o 

Richard E. Turner, Thomas H. Rees, 
and Gerard E. Woodbury 
Langley Research Center 

SUMMARY 

A new primitive -variable general circulation model of the ocean has been developed 
in which fast external gravity waves may be suppressed with rigid-lid surface constraint 
pressures. The surface pressure filtering technique is easier to apply than the conven- 
tional stream function technique, and since the resulting model is derived in terms of 
primitive variables, it is conceptually simpler. From theoretical considerations, the 
model appears to be a candidate for limited region simulation, as well as simulation of 
global regions. 


INTRODUCTION 

Ocean circulation models are generally formulated either in terms of a stream 
function, such as the Bryan model given in reference 1, or else as a primitive -variable 
model, such as the one developed in reference 2. Both types of models have problems 
that one might like to avoid. 

Stream function models are difficult to apply to multiconnected regions involving 
islands. Since the momentum equations are formulated in terms of a stream function, 
such models may be conceptually complex; and since stream function models simply can- 
not represent divergent flow fields, they are limited to large-scale problems where tidal 
motions can be neglected. (Tidal motions become more important in coastal zones and 
tidal basins.) Primitive -variable models, on the other hand, are conceptually simple and 
flexible. Islands and divergent flow fields are easily represented with primitive -variable 
ocean models. Unfortunately, primitive -variable models which do not filter fast external 
gravity waves associated with free -surf ace topography require large amounts of computer 
time. 

The present model is a pseudo-primitive-variable model in which the fast external 
gravity waves are filtered by applying constraint pressures computed from the rigid-lid 
assumption. The resulting model can use a long time step and yet has the conceptual 
simplicity normally associated with a primitive -variable model. The model's vertical 



1 1 III II 


structure is represented by a generally scaled variable that allows a simple but flexible 
treatment of bottom topography and land-sea boundaries as well as free -surface topog- 
raphy in situations where free-surface geometry is important. Islands are easily repre- 
sented in the model because the boundary conditions at the islands are known a priori as 
opposed to stream function models wherein boundary values of the stream function around 
islands are difficult to obtain. 

The surface pressure technique for filtering the external gravity wave is not entirely 
new. A similar approach is presented in reference 3; however, the present model has 
several advantages. First, the mathematical model presented herein allows a general 
stretching of the vertical coordinate rather than the linear stretching given in reference 3. 
Second, the model is derived in terms of spherical polar coordinates and is suitable for 
large-scale circulation problems. The rigid-lid model has been investigated analytically 
for open boundary situations and appears to be a candidate for such problems. Finally, 
the subgrid mixing formulation is presented in terms of a strain-rate tensor that properly 
vanishes when the fluid rotates in the rigid body mode. Consequently, erroneous shear 
stresses are not produced by simple rigid body rotations as in most previous works. 

Thus, the present model fulfills a basic need that exists in the field of dynamic 
oceanography. A conceptually simple model has been established that can operate effi- 
ciently with a long time step. It is well suited to large-scale closed region simulation 
where free-surface topography is not important, as well as to small-scale open lateral 
boundary simulations where free-surface topography may be important. 

NASA is presently studying potential use of satellites to monitor water pollution. 
Efforts will be made to determine the enhancement of remotely sensed data by use of 
circulation models. Since there is no way at the present time to penetrate the ocean 
depths by remote sensing, the initial attempt to study ocean pollution problems by remote 
sensing must be done by inferring subsurface phenomena from circulation models used in 
conjunction with remote and in situ measurements. 

A description of the mathematical basis for the present model is given in refer- 
ence 4. In the present report, the specialization of the mathematical model for ocean 
simulation is presented as well as preliminary computations on a 10° global grid. A sim- 
plified analysis of boundary condition specification for computation of rigid-lid constraint 
pressures is given in appendix A by Huw C. Davies. Detailed calculation for a model of 
the surface layers of the North Atlantic Ocean is presented and discussed in reference 5. 
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SYMBOLS 


A 

B 

C D 

c v 

D(x“,t) 


D f 

E 


f 


g ii 

H(xi,t) 


I,J,K 

K H 

Kq 

K v 

L iJ 


determinant of quasi-horizontal metric tensor in spherical polar coordinates; 
also, area in appendix A 

parameter in equation of state 

drag coefficient for bottom friction calculations 

specific heat at constant volume of sea water 

array specifying total depth distribution at horizontal position (x\x^) 
and time t 

flat -top depth distribution array or depth relative to the mean sea level geoid 

internal energy per unit mass 

group of terms in equations (24) and (25) 

Coriolis parameter 

acceleration due to gravity 

covariant element of the metric tensor 

four -dimensional function specifying instantaneous height of horizontal para- 
metric surfaces referred to the mean sea level geoid 

grid point indices 

horizontal subgrid mixing coefficient 
Von Karman constant for turbulent flow, 0.4 
vertical subgrid mixing coefficient 

physical components of the deviatoric strain-rate tensor 
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cm 


1 square root of second scalar for f 1 . 

J 

a i 

£ j mixed tensor components of deviatoric strain-rate tensor 

normal unit vector 
P pressure 

P approximate pressure 

P £ surface constraint pressure 

Pf hydrostatic pressure for rigid lid 

q perturbation constraint pressure in appendix A 

r mean radius of Earth 

S salinity mass fraction 



S J 


s 

T 

t 


u i 


U 


a 


U 


"x3=0 


V J 


source term 

physical velocity component of coordinate grid point in jth direction relative 
to rotating Earth 

arc length 

temperature, °C 

time 

physical component of fluid velocity vector relative to rotating Earth 
horizontal velocity component at end of first stage of integration 
horizontal velocity at ocean bottom 

physical component of fluid velocity relative to coordinate grid points in 
jth direction 
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X 3 


e 

X 


v 


4 

P 



T 


$ 

U> 

CO 

Indices: 

U 

t 


total depth of ocean in transformed coordinates 
contravariant coordinate variables of reference coordinate system 
Cartesian coordinates 
grid increment 

horizontal perturbation velocity in appendix A 
generalized density 
parameter in equation of state 
elemental fluid volume 

function specifying layer thickness distribution 
density of sea water 

density of sea water at standard conditions 

elements of two-dimensional square matrix used to compute Coriolis force 

time increment 

general physical parameter 

angular velocity of Earth 

mean molecular weight of sea water 

take on values 1, 2, or 3 
computed at time t 
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take on values 1 or 2 only 


ot,l 3 

1,2,3 particular coordinate direction 

The summation convention is used for multiple indices unless otherwise indicated. 
Notation: 

( y space and time average 

( ) finite increment in ( ) 

Primes denote the difference between instantaneous value and space-time average. 

ANALYSIS 

The model is based upon conservation of mass, momentum, salinity, and internal 
energy. The Navier-Stokes equations for a rotating, nearly spherical Earth are somewhat 
simplified by several important assumptions and approximations. First, the equation of 
motion for the vertical velocity is reduced to the hydrostatic approximation by neglecting 
local acceleration and other terms of the same magnitude. Second, the Boussinesq approx- 
imation, which neglects density differences except in the buoyancy term, is made. Third, 
since the equations are solved numerically on a finite -difference grid, stresses and pro- 
cesses occurring on scales too small to be resolved on the grid are parametrized or 
neglected. Molecular viscosity and conductivity are neglected. An empirical equation of 
state, giving density as a function of temperature, salinity, and pressure, is adopted. 
Finally, the "rigid-lid" approximation, which suppresses external inertia-gravitational 
oscillations by forcing the free surface to conform approximately to the mean sea level 
geoid, is optionally incorporated. 


Coordinate System 

For flexibility, the equations are formulated in spherical polar tensor notation. The 
coordinate system chosen for the ocean model is described in detail in reference 4. The 
coordinates are divided into two horizontal components x and x and one vertical 
component x^. The bottom horizontal parametric surface (x^ = o) follows the ocean 
bottom topography, while the upper coordinate surface (x^ = X^) follows the upper surface 
of the ocean. Intermediate surfaces can be located as desired. To add flexibility to the 
model vertical structure, the equations of motion are formulated in terms of density mul- 
tiplied by the vertical scale factor to form a "generalized density" 
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( 1 ) 


^ “ p o^ g 33 



‘) 


The family of x 3 -surfaces is controlled by the choice of which can vary with x*, 

x^, x^, and time. An x 3 -surface is constructed by translating from x 3 = 0 a physical 

px 3 

distance j dx perpendicular to the local x J -surfaces. 


The x* parametric lines at the surface x 3 = 0 are intersections of constant- 
longitude planes with the ocean bottom. The x^ parametric lines at the surface x 3 = 0 
are intersections of constant-latitude cones with the ocean bottom. Since the "horizontal” 
velocity components follow the x A and x parametric lines, these components are gen- 
erally not perpendicular to the local gravity vector. Therefore, in the momentum equa- 
tions the components of gravity along the quasi -horizontal coordinate lines are included in 
the momentum balance. The validity of this formulation presupposes a small angle 
between dx 3 and the local gravity vector. 


From reference 4, the square of the elemental arc length is shown to be 
approximately 


(ds Y 


» (r dx*) + (r sin x* dx^) + ggg (dx 3 ) 


( 2 ) 


so that 


*11 * W' 


g 22 


~ (r sin x*) 


( 3 ) 


The determinant of the quasi-horizontal metric tensor is 
A = g ll g 22 

12 Q 

The unit vectors along dx , dx , and dx° form an orthogonal triad, and since ggg is 
allowed to vary with time, the grid points have physical velocities Sp Sg, and Sg rela- 
tive to rotating Earth. The momentum equations are presented in terms of physical veloc- 
ities Up Ug, and Ug relative to rotating Earth such that 
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where is the physical velocity of the fluid relative to the grid. All velocity com- 
ponents are given with respect to the local (dx^,dx^,dx^) triad. In the model presented 
herein, and S 2 are always neglected as in reference 4, and S 3 is determined by 
the manner in which the vertical structure is chosen. The technique used to compute V 3 
is presented in a subsequent section of this report. 

Dynamic Equations 

A rigorous derivation of the governing equations in spherical polar tensors is pre- 
sented in reference 4. The equations are given here in final form only. The continuity 
equation in terms of the generalized density £ is 


at 


_1 a_ 

y/A dx a 





(4) 


where the second and third terms represent horizontal and vertical divergence, 
respectively. 

The dynamic equation for momentum conservation in the -direction (north-south 
direction) is 


_a_ 

at 




./l33 sp^ 2fu „ CQS x _ 

^n 8xl 2 


,ys 77 ax 1 /a l 


<VV>' 

'J^aa 


9x ' v^33 / \ \[%n \fiii 


r 9Sq 

S TT & 


■ u. 


JzTi 39x1 


(5) 


In the above equation, the three terms on the left side (from left to right) represent 
accumulation and horizontal and vertical advection. The first term on the right side is 
the pressure gradient term. The second term on the right side is the Coriolis force term, 
where c 0 is the angular velocity of the Earth. The third term is the oblique gravity term 
which arises from departure of the quasi-horizontal parametric surfaces from the true 
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local horizontal. The instantaneous height of the horizontal coordinate surface referenced 
to the mean sea level geoid is 



where Dj is the depth relative to the mean sea level geoid. The fourth and fifth terms 
on the right side of equation (5) represent horizontal and vertical subgrid diffusion, respec- 
tively. The brackets ( ) indicate space and time averaging, and the primed quantities 
indicate the difference between the instantaneous value and the space-time average. The 
next to last term represents centrifugal force, while the last term arises from motion of 
the vertical grid structure. Similarly, the momentum equation in the x^ -direction (east- 
west direction) is 



The momentum equation in the x^ -direction is discarded in favor of the hydrostatic 
approximation. 

The equation for conservation of salinity is 




a 


actj 





a ( W)\ 

9x3 \ ^33 ) 


(7) 
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where the three terms on the left side represent accumulation and horizontal and vertical 
advection, and the second and third terms on the right side represent horizontal and ver- 
tical subgrid transport. The term S A %s) is a source term. Similarly, the equa- 
tion for conservation of internal energy E is 



where work of compression is ignored. 


Density and Pressure Computations 

Since the density of sea water is numerically very close to 1 g/cm 3, the Boussinesq 
approximation which neglects variations in p, except where multiplied by g, is a rea- 
sonable one and is used herein. This approximation simplifies the calculations. For 
example, the generalized density £ is given simply by £ = ^ ggg. However, since the 
actual density is required in computing pressure, the technique for computing p is out- 
lined in this section. 

The empirical equation of state from reference 1 which gives density p (in g/cm^) 
as a function of sea pressure P (in bars (1 bar = 0.1 MN/m^)), salinity mass fraction S 
(in ppt), and temperature T (in °C) is 


p( S,T,P) 


B 

1.000027(A + 0.698B) 


( 9 ) 


where 


B = P + 1 + 5890 + 38T - 0.375T 2 + 3S 


and 


A = 1779.5 + 11.25T - 0.0745T 2 - (3.8 + 0.01T)S 
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Temperature T is given by 


T = — 
c v 

where C y is specific heat at constant volume and u) is mean molecular weight of sea 
water. The seas are assumed to be in hydrostatic equilibrium, so P is computed by 
integrating (with respect to x^) the hydrostatic equation, 

- - g e#33 (10 > 

For given S and T, equations (9) and (10) form two simultaneous equations in two 
unknowns. To simplify the solution, an approximate pressure P is first computed from 
equation (10) by setting p = 1 g/cm 3 . The density is then taken as p = p(S,T,P) in 
equation (9), and equation (10) is again applied to give a closer approximation to P. 

Vertical Velocity 

In equations (4) to (8), Vg is needed to compute vertical transport. The formula- 
tion allows ^ ggg to be chosen quite generally provided that the following boundary con- 
ditions at the bottom and upper surfaces are met: 

Vg = 0 at x 3 = 0 and x 3 = X 3 (11) 

In the present model, Vg (and hence S 3 ) is determined by requiring that 

? = E>(x a ,t) |(x 3 )p q (12) 

where the function D specifies the total depth distribution. The function i;(x 3 ) which 

O 

defines the x profile of the generalized density £ is restricted so that 



Otherwise, |(x 3 ) can be any positive, single-valued function. Once this function is 
chosen, the expression for Vg can be found by substituting equation (12) into equation (4), 

integrating with respect to x 3 from the bottom surface (x 3 = 0 ) to some specified 
x 3 -surface, and applying the boundary conditions of equation (11). The result is 
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( 14 ) 



An expression for 9D/9t can be obtained by evaluating equation (14) at x 3 = X 3 where 
V 3 = 0 (eq. (11)) and solving for 9D/9t; that is, 



Equation (14) results directly from specifying £ by equation (12). The determination 
of V 3 in this manner allows one to specify the x 3 grid point as a given fraction of the 
total depth regardless of the spatial variation of the depth. It should be pointed out that 
because of the formulation of the governing equations in terms of moving grid points, 
other methods of specifying Vg could be adopted. For example, a hybrid Euler -Lagrange 
system of coordinates could be achieved by specifying Vg = 0 along with S.^ = Sg = 0. 

Gravitational Instability 

In the ocean, convective overturning maintains a stable vertical density profile. The 
hydrostatic approximation employed herein neglects vertical accelerations and prevents 
the model from simulating convective overturning. Therefore, a stable density profile 
must be maintained in the model artificially. Every water column is tested separately 
(at each integration step) for a gravitational instability. First, the reference-pressure 
densities are computed for each grid point in a given vertical column. If the grid values 
of reference -pressure densities decrease monotonically with grid distance from the 
bottom, then no gravitational instability exists in the given column, A gravitational insta- 
bility is detected at level K (where K denotes the number of grid points from the 

bottom) when p p (K) = p p (K+l). The instability is alleviated by mixing 

^surface ■'surface 

the two nodal cells to uniformity in temperature, salinity, and velocity. The new values 
of temperature, salinity, and velocity for both cells are computed to be 

™ - TfIC4.11 - g(K) T ( K > + ^ K+1 > T ( K+1 ) 

T(K) - T(K+1) K) - (K+1) 


S(K) = S(K+1) = 


5(K) S(K) + ?(K+1) S(K+1) 

t(K) + 5(K+1) 
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Ill 


U a (K) = U a (K+l) = 


?(K) U a ( K) + £(K+1) U g (K+l) 
?(K) + ?(K+1) 


The total internal energy, salt, and momentum in the two cells are, therefore, conserved. 
The readjustment scheme given is chosen for simplicity and is intended to serve as an 
example only. A more complex and rigorous technique is given in reference 1. 


Subgrid Mixing 

The action of the unresolved subgrid waves on the resolvable waves is treated by a 
nonlinear kinematic eddy viscosity in the horizontal directions and by a specified viscosity 
in the vertical direction. 

The subgrid velocity correlation tensor is taken to be proportional to the strain-rate 
tensor; that is, the velocity correlation components are related to the tensor strain-rate 
components by (from ref. 4) 

< U c/V) = " K H^Q!i3 < U o/ U 3’> = ~ K V^O!3 


where 


Lij = fi 


j( g ii/ g jj) 


1/2 


(i and j not summed) (16) 


and the t 1 ^ are given for the quasi-two-dimensional case as 


_9 

ax 


i ( u i sin xl ) 


au 


8x 


Y - 2 cos Uj 


// g 22 


1^2 ~ s ^ n 


9Ui a / i\ i 

— ^ + — -IU 0 sin x 1 ) - 2 cos x 1 U 9 
8x 2 9x 1W ' 


A/ g 22 


(17a) 


(17b) 


.1 9U X 

£ o = 


3-^37 V^l 


(17c) 


£ 2 j = tJ(sm x 1 ) 


(17d) 
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(17e) 


-2 -1 
1 2 = -*■ 1 


£ 3 = 


8x3//^ 

(17f) 

^3^11/^33) 

(17g) 

^ 2 3( g 22/ g 33) 

(17h) 

0 

(171) 


In the horizontal directions the strain-rate tensor is the basis for subgrid mixing. The 
strain-rate tensor, computed by assuming quasi-two -dimensional motion (zero strain rate 
in the vertical direction), is used to compute the horizontal kinematic eddy viscosity Kjj. 
The tensor analogy for K^, modeled after reference 6, is from reference 4, 

K H =/2(K 0 A ) 2 f (18) 

where K 0 = 0.4, A is the reference grid dimension, and 


£ = 



(19) 


In the vertical direction the subgrid mixing coefficient Ky is specified and may vary 
with depth. 

The horizontal subgrid mixing for salinity and sensible heat flux from reference 4 
are 


< u «'s> 

<V E ') = 





K H C V 3T 



(o 

not summed) 

(20) 

(a 

not summed) 

(21) 
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where C v is the specific heat at constant volume and a; is the mean molecular weight 
of sea water. The corresponding vertical subgrid mixing for salinity and sensible heat 
flux are 


< U 3' S ’> - 


Kv S/v g33 

(22) 

K V C v _9T_ 
/S33 " 9x3 

(23) 


The approximation used herein is that Ky is identical for velocity, salinity, and temper- 
ature. If one desires, different values of Kjj can be used. 

Surface Pressure Approach to Rigid-Lid Approximation 

For the present model a simple technique has been developed to filter external 
gravity waves. Incorporation of this filtering technique allows the model to simulate 
large-scale flow efficiently without troublesome external gravity waves. However, the 
model formulation allows deletion of the external gravity wave filter to simulate tidal 
motions and free-surface topography effects in small-scale studies. The gravity wave 
filter discussed in this section (which is a form of the rigid-lid approximation customarily 
used in stream function models of large-scale ocean circulation) is applicable both for 
problems having closed natural boundaries and for problems having open boundaries. 

The technique is based upon the idea that fast external gravity waves can be filtered 
from the model by computing and applying the pressure that a rigid lid would exert on the 
ocean surface to negate time -dependent depth fluctuations. Application of the constraint 
pressure involves only the horizontal momentum equations. 

The constraint pressure P c is introduced into the momentum tendencies given in 
equations (5) and (6) by replacing P with P^ + P c , where P f is the hydrostatic pressure 
based upon the flat-top depth. By grouping terms, equations (5) and (6) can be rewritten as 


_d_ 

dt 




\/ g 33 9P c 

fiu 9x1 


V^33 9P c 

\f & 22 9x2 


+ / g 33 F l 
+ /S33 F 2 


(24) 


(25) 
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where F^ and F 2 contain all the remaining terms ^divided by ^ggg j in equations (5) 

and (6). Equations (24) and (25) contain three unknowns, ^^ U 2^’ arid P c' 

The third equation needed to close the system can be obtained after noticing that P c 
must be consistent with the requirement that the depth not vary with time. From equa- 
tions (1), (12), and (13), the ocean depth given at a point (xl,x 3 ) is 



dx 3 


(26) 


Thus, P c must be consistent with 


P 


3D 

o at 



at 


dx 3 = 0 


(27) 


Substituting equation (4) into equation (27), integrating, and applying the boundary condi- 
tions on V 3 (eq. (11)) result in 


P 


3D 

o at 


_1 a_ 

\JX dx a 





dx 3 = 0 


(28) 


When equations (24) and (25) are substituted into the derivative of equation (28) with 
respect to time and equation (26) is used, there follows 


1 a y/A 



dx 3 = + 


_i a ( \Kd 9 P c\ 

y/A ax“ \ ^ aa 9x a / 


_1 a_ 

s/A dx a 





F 


a 


dx 3 


which because of equation (27) vanishes leaving 


J a_. 

A/ad 9 P c\ 

sJa ax a 



_1 a 'JX 




(29) 


as the governing Poisson equation for P c> Equation (29) is correct in theory; however, 
in practice P cannot be computed exactly and additional requirements are needed to 
insure that 8D/at vanishes. These additional considerations are given in a subsequent 
section which deals with the numerical details. 
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Boundary Conditions 


Boundary conditions to the dynamic equations are required at the ocean floor, along 
land-sea boundaries (and open boundaries, if any), and at the air-sea interface. Those 
boundary conditions which have not already been discussed are noted in this section. 

Ocean bottom friction.- For currents in the bottom layer, friction terms are com- 
puted by use of a neutrally stable Prandtl layer model, originally used in an atmospheric 
circulation model (ref. 6). In essence, the Prandtl layer model assumes that vertical 
subgrid transport is in equilibrium within a thin bottom layer. The momentum transport 
in the Prandtl layer is derived from reference 6 as 



where 


U 


a 


= 0 and 


x 3 =0 


Cp is a drag coefficient. 


(30) 


Air-sea interface.- The interaction terms at the air-sea interface are wind stress 

p ^U^'Ug'^, heat transport rate p ^E'Ug' and salinity counterflux p^S'Ug'^. The 

effect of evaporation upon the water mass of the system is neglected. These interaction 
terms must be specified. Methods for computing them are given in reference 6. 

Land-sea boundaries.- The boundary condition at coastlines requires that mass 
transport normal to the beach be zero at the boundary. This condition is automatically 
satisfied by requiring the depth to be zero at boundary points. Thus, the coastline always 
follows a parametric line connecting boundary points. 

When the rigid-lid approximation is used, an additional boundary condition - that 
the mass transport normal to the beach in a water column adjacent to a land boundary is 
zero - is required. This condition is imposed by subtracting the mean values of Uq, 
throughout the column. 

Open boundary conditions for rigid-lid applications. - A simplified analysis of bound- 
ary condition specification for the momentum equations and surface constraint pressure 
is given in appendix A. The approach is a version of a method by Davies in reference 7 
for studying the uniqueness of solutions to the shallow water equations in a finite region. 
Along inflow boundaries, the velocity components U a must be specified. The surface 
constraint pressure can be computed from an equation similar to equation (Bll) in 
appendix B. Along outflow boundaries, either normal outflow velocity or else surface 
constraint pressure must be specified. When velocities are specified, they must satisfy 
mass conservation. Since P c has an analogy with the surface topography caused by 
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external gravity waves, one might specify P c by measuring the free-surface topography 
along the outflow boundaries. In many oceanographic flow situations, the geostrophic 
approximation in which the Coriolis force is balanced by the pressure gradient is justified. 
In the case of a limited region, it should be possible to choose boundaries through regions 
where the geostrophic approximation is valid. Then along an outflow boundary, approxi- 
mate values for P c would be computed from outflow velocities. To see this, recall 
equation (5) and retain only the pressure gradient, the oblique gravity term, and the 
Coriolis term. The sum of all three terms vanishes in the geostrophic approximation; 
thus, if dx* lies along the outflow boundary, one has the expression 

- -^51 + 2 £o>U 2 cos x 1 ~ 0 (31) 

^IT 8x \l%n 8x 

which can be integrated along the outflow boundaries to specify P c . 

COMPUTATIONAL SCHEME 

The dynamic equations are solved over a three-dimensional grid by a finite- 
difference scheme which is similar to the method used in reference 6. The technique, 
central differencing in time and space, is often referred to as leapfrog. 

Lattice Structure 

The dependent variables are laid out over the three-dimensional lattice so that 

^(xW 3 ) = <I>(lAx 1 ,JAx 2 ,KAx 3 ) = $(I,J,K) 

where the Ax A translation is from north to south, Ax* from west to east, and Ax 
is along the local x^ parametric line, perpendicular to the ocean bottom and positive up. 
On the two-dimensional horizontal grid, all dependent variables are defined at time 
t - r on half of the grid points and at time t on the remaining grid points. The time 
structure does not change with x 3 ; in other words, at latitude index I and longitude 
index J, the variables are defined at the same time for each point in the column. Thus, 
the twin time lattices may be referred to as "even" or "odd,” depending upon whether 
I + J is even or odd. The even lattice, defined at time t, for example, is used to update 
the odd lattice from time t - r to t + t. Then the odd lattice at t + r is used to 
update the even lattice from time t to t + 2 r. 
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Space Differencing of Dynamic Equations 

To illustrate the computational sequence, integration of the salinity conservation 
equation is described. For simplicity put 


^ “ \l s aa~ 33 ~ ? " 1 


(v/a ?s) = o 


in equation (7) and get the simpler result 


as _ _ _a 

3t ax 


H( SU a )-^( SV 3)-^< S ' U -')-i3< S,U 3'> 


The subgrid mixing given by equations (20) and (22) with the assumed simplifications 
becomes 


< s 'V> = - k h 


< S ’ U 3’)' K v5 


The quantities Vg and ^S'Ug'^ are computed at (I,J,K±l/2), whereas the quantities 

<S'U a ’) are computed at (Itl/2,J,K) and (I,J±1/2,K). Thus, to update dependent 

variables on the odd lattice from time t - r to time t + t, one uses the following 
procedure: 

(1) Use the even lattice at time t to compute 

a . S(I+1,J,K) ^(1+1, J,K) - S(I-1, J,K) U 1 (I-1,J,K) 

dx a v 01 ' 2 Ax* 

S(I,J+1,K) U 2 (I,J+1,K) - S(I,J-1,K) U 2 (I,J-1,K) 

2 Ax 2 


i( SV 3) 


S(I,J,K+l/2) V 3 (I,J,K+l/2) - S(I, J,K-l/2) Vg(I,J,K-l/2) 


ill! I Ill I 


I I I I II II Mill 


III I III II I 


II 


I I II III I I I I II I 


where 


S(I J K±l/2 ) = S(I-fl,J,K±l) + S(I-l,J,K:bl) + S(I ? J-hl^ K±l) + S(I,J-1 ? K±1) 

7 ? g 

, S(I+1,J.K) + S(I-1,J,K) + S(I,J+1,K) + S(I,J-1,K) 

8 


(2) Use the odd lattice at time t - r to compute 



Ky(K+l/2) 

S(I,J,K+1) - S(I,J,K)“ 

Ky(K-l/2) 

j V 

|S(I, J,K) - S(I,J,K-1)1 

Ax 3 

_ Ax 3 

Ax 3 

[ Ax 3 J 


Recall that Ky is specified at each vertical level. 

(3) Use both the odd lattice at time t - r and the even lattice at time t to compute 


dX 01 ^ Ua ) ' dx a ( KH dx a 


S K H (I+1 / 2?J;K) [ S(I+1,J,K) - S(I,J,K)1 + K H (1 ' 1/2?J>K) f s(I,J,K) - S(I-l,J,K) j 

Ax 1 L Ax 1 J Ax 1 L Ax 1 J 

_ Kn^J+lAK)^ J+1 K) _ S ( lfJ>K n + K h (I,J- 1/2,K)^ (IiJ K ) _ g ( It j_ ltK ) j 

Ax 2 L Ax 2 J Ax 2 L Ax 2 J 


where, for example, K H (I+1/2,J,K) is computed from equations (16) to (19). From equa- 
tions (16) and (17), let 



where 


U 2 (I+1,J,K) + U 2 (I+1,J+1,K) + U 2 (I,J+1,K) + u 2 (i,j,k) 

A ( U 2/ = 4 

and 

U 9 (I+1,J,K) + U„(I+1,J-1,K) + U«(I,J-1,K) + U«(I,J,K) 

b (U 2 ) - -A 2 _ ? 2 

Then 

. . U 2 (I+1,J,K) - U 2 (I,J,K) A(Uj) - Bp!) 

L 12 = L 21 = + ^2 


where A(UjJ and B^UjJ are defined in a manner similar to A(Ug) and B(Ug), 
respectively. Using equations (18) and (19), one gets 

f~ 2 

K h (I+1/2,J,K) = 2(Kq) 2 Ax* Ax 2 (L n ) + (L 12 
The values of Kjj(I-1/2,J,K), Kjj(I,J+ 1/2,K), and Kjj(I,J- 1/2,K) are computed by 

A A 

simple permutations of the indices. Note that and L^ 2 were simplified because 

of the assumption y/K = ^ Jg ~ = ^g 33 = 1. 

Time Differencing for Dynamic Equations and Constraint Pressure 

Since the constraint pressure affects only the momentum tendencies, it is convenient 
to break the integration process into two distinct stages. In the first stage the dependent 
variables £, S, and E are updated with equations (4), (7), and (8), respectively. Since 
P c is not known in the first stage, and U 2 are updated with P c = 0 in equa- 

tions (24) and (25). At the end of the first-stage integration, £, S, and E have the cor- 
rect new values while Uj and U 2 have incorrect values denoted by and Jg. For 

clarity, equations (24) and (25) may now be rewritten in terms of Uj and Ug by using 




21 


( 35 ) 


(?u 2 ) l+T (CU 2 )‘' T , (5U 2 ) t+T - (c t+T u 2 ) ( (c t+T u 2 ) - (CU 2 ) t - T 

2t " 2r 2 r 2r 


The first-stage integration of the momentum equations is defined to be 


(? t+T Uj) - (?u 1 ) t_T _ 


2t 




(36) 


and 




2 T 


\/^33 ] 


(37) 


In the second stage of integration, the correct values of and U 2 are found from the 

equations 


(?u x ) t+T - 



2t 



and 


(gU 2 ) t+T - (^ t+T U 2 ) _ \^3P C 

2r V^22 9x2 


(38) 


(39) 


In order to solve equations (38) and (39) for Uj^ +T and U^ 4 " 7 , a third equation 
for P c is needed. In theory, one can substitute equations (38) and (39) into the time 
derivative of equation (28) and solve a Poisson-type equation for P c at time t + r as 

denoted in equation (29). Then Uj^ +T and U 2 t+T can be found from equations (38) 
and (39). The fallacy in this simplistic approach is that P c cannot be found exactly; 
consequently, equation (27) would be violated over a long period of time. The proper 
approach, then, is to apply equation (28) in such a way that the time variation of D(x a ,t) 
is bounded. More specifically, equation (28) is replaced by its time -discretized counter- 
part evaluated at time t + t, so that 
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(p 0 - 9 < 40 > 



l a_ 

sfA dx a 




dx'- 


Since P„ cannot be computed exactly, D t+2T and D* must differ from the rigid-lid 
depth D^. In the order of computation, D will have been determined by previous com- 
putations, but the value of T^+It can still be specified. In practice, D t+3T is set 
equal to in equation (40). Substituting equations (38) and (39) into equation (40) yields 


D f -D' + x 

1 

a 

Sa 

(2t) 2 2t 


dx a 

\J^aa 

w 



dx 3 


1 a I 

/ADf ap\ 

V'A dx a 

l g a« dx a 


(41) 


One can obtain equations for the surface pressure along open inflow and closed lateral 
boundaries from equation (41). This is done by applying the fictitious boundary conditions 
given in appendix B as equation (B6) for closed boundaries or equations (B12) and (B13) 
for inflow boundaries. Equation (41) can be solved for P c by relaxation when the bound- 
ary conditions for P £ are known. 

In the derivation of equation (41), no account was taken of the possibility that F a 
in equations (24) and (25) might contain U c J i+T , as would occur when the Coriolis term is 
integrated implicitly. Consequently, equation (41) can be solved in a straightforward 
manner only when an explicit integration technique is used. In ocean circulation problems, 
when long time steps are achieved by filtering external gravity waves, the minimum upper 
bound for the time increment is determined by the Coriolis term (fr SI). Since implicit 
treatment of the Coriolis term increases the allowable time step, a simplistic technique 
has been devised for implicit integration of the Coriolis term in conjunction with the sur- 
face pressure constraint technique for steady-state solutions. Equations (24) and (25) are 
rewritten to include values of P c from the previous step, so that 

t_2 T 

JLfqj \ = a ( AP e) + /g— F - tEs. ia not summed) (42) 

att a) yftZ 9X “ V ^ 9X “ 


f 4-0 t-r 

where AP C = P c - P c . The first stage of the integration is then defined by the 
relationship 


(C t+T D a ) - (CD a )‘- T 

2 t V g 33*of 


g 33 9P c 


t-2r 




act 


ax 1 


a 


(a not summed) (43) 
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and with the Coriolis force computed as 


The second integration stage is performed with 

(cu a ) t+T - (; ttT u tt ) #^ 8 ( ap c) ( 

2t 9x “ 

When U a t+T converges to a steady-state solution, P c t converges to a steady-state 
value while converges to U Q , t+T - The Coriolis force is thus treated implicitly. 

MODEL SIMULATION RESULTS 


A computer model based upon the mathematical model described herein was 
developed on a global 10° grid, and some preliminary numerical simulations have been 
performed. Since a 10° resolution is too coarse to be compared in detail with the real 
world, only qualitative features of the flow field are discussed. 

In particular, a one -layer flat -bottom model was developed with continental geometry 
and surface wind stress included; salinity and heat transport were not included because of 
the long time scales required for these quantities to reach quasi-equilibrium conditions. 
The annual mean wind stress used to drive the flow field is plotted in figure 1 and was 
obtained from reference 8. The model was started from rest by applying the wind stress 
and calculations continued until quasi-equilibrium was reached. The resulting flow field 
is presented in figure 2 with the zero meridian at the left boundary and extending to 360° 
at the right boundary. The model simulation area runs from the North Pole at the upper 
boundary to the South Pole at the lower boundary. The real-world flow field, taken from 
reference 9, is shown in figure 3. The model flow field contains large gyres in the cor- 
rect positions and rotating in the proper sense but without realistic detail. The equatorial 
currents are observed to be weak. The western intensifications are not present except in 
the North Atlantic, and the west wind drift is observed in the South Pacific. The probable 
cause for loss of realism in the model flow field is the coarse grid spacing and the single 
layer for vertical resolution. During the computer experiment, the ocean's depth remained 
constant to five significant figures; therefore the constraint pressure technique was 
working properly. From this experiment, it was concluded that the overall approach is 
valid and ready for further development on a finer scale. 
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For detailed simulation results, the reader is referred to reference 5 wherein the 
results of a model of the North Atlantic Ocean on a 2.5° grid for wind -driven surface 
currents are presented and compared with real-world data. 

CONCLUDING REMARKS 

A new pseudo-primitive -variable ocean circulation model has been developed along 
with a surface pressure constraint technique for filtering fast external gravity waves. 

The coordinate system follows the ocean bottom and surface. The advantages of the 
mathematical model presented herein are (1) the ease of determination of boundary con- 
ditions, (2) the ability to either resolve or filter external gravity waves, and (3) the ability 
to utilize coastal bottom topography. The model has been tested under some simple situ- 
ations which indicate that the surface pressure constraint method can produce divergence- 
free flow fields in circulation problems involving wind stress. The results of the present 
report indicate that the new ocean model warrants further investigation. 

Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
July 15, 1976 
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SPECIFICATION OF BOUNDARY CONDITIONS FOR COMPUTATION 
OF RIGID-LID CONSTRAINT PRESSURE 


Huw C. Davies 
University of Reading 
Berkshire, England 


In this appendix an idealized analysis is made of the specification of boundary con- 
ditions for the constraint pressure field. The approach is similar to, but simpler than, 
the technique used by Davies in reference 7 to specify boundary conditions for the shallow 
water equations. 

The analysis begins with the equation of motion for a constant-density fluid confined 
between two parallel plates with a uniform separation distance and bounded by inflow, out- 
flow, and no-flow surfaces. The constraint pressure gradient and the Coriolis force are 
taken to be the sole driving forces; thus, equations (5) and (6) simplify in rectangular 
Cartesian coordinates to 


au 


a 


at 






)■ 


8P, 


9z 


t , +1 °a.A 


(Al) 


where 



1 

0 


The method of solution is to find the rate of growth for perturbations of U^ and then to 
specify P £ so that perturbations in U^ are minimized. Let perturbations in U^ 
and P be and q, respectively; thus, equation (Al) becomes 

^ Ct 


^( U a +e J + ^|(V e /3)( 


Wa + e a 


) " ( Pc + 


C l) +fa « /3 ( U / 3 +6 / 3) < A2 > 


Upon neglecting quadratic perturbation terms and using equation (Al), the governing equa- 
tion for e„ becomes 

a 
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8e 


a 

at 




_ 3q 


9z 


a +f W/3 


(A3) 


Next, a specification of q is sought which causes the volume integral of 6 ff e ff to 
decrease in time and thus minimize the total perturbation energy. An equation for e a e a 
is found by multiplying equation (A3) by e a to get 


_a_ 

at 




+ e / 3 U 



-€ 


a 



(A4) 


The Coriolis force cancels in the summation. Then the time derivative of is 

integrated over a volume v to yield 


d_C / e a 6 a\ 
9t J if \ 2 / 


d„=.r e -i- 

/U 0 e ) 

\ a dzP 

l 0 a ) 

- f ^ aq 

dv 

01 dz a 



dv 



(A5) 


In order to analyze the effects of boundary conditions on the left side of equation (A5) ? the 
first and last terms on the right side can be rearranged to get 


_a_ 

at 



dp = - 



V«'« 


) 



dz^ 



f — ^ — (e ^q') dv + f q — — dv 
^vdz a ' a ’ Jii 


v dz L 


(A6) 


Consider the right side of equation (A6), wherein 

= o 

az^ dz@ 

The first two integrals combine to give 
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the third integral simplifies to 

and the last term vanishes leaving equation (A6) as 



Next, suppose that initially 6 ff e a vanishes everywhere on the interior of v. Then only 
two surface integrals remain on the right side of equation (A7), those being 


- Ji K and -J A M ,, <. dA 


where n^ is positive outward on surface dA of v. Consequently, equation (A 7) 
simplifies to 



where dAp dA Q , and dA c refer to inflow, outflow, and closed boundaries, respec- 
tively. Since the first term on the right side of equation (A8) increases the perturbation 
energy, e e /2 = 0 must be maintained on inflow boundaries to insure a stable com- 
putation. The second term tends to decrease perturbation energy so that /2 may 

be nonzero on outflow boundaries. Since e Q; e Q ,y^2 must vanish on inflow boundaries, the 
third term vanishes, as well as on closed boundaries. Thus, equation (A8) 

reduces to 
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A f e a e « 

at J„ 2 


dv 


*-y A 

A o 


dA 0 


(A9) 


The vanishing of on inflow boundaries is equivalent to specifying on 

inflow boundaries. Because of equation (A9), 6 a vanishes along outflow boundaries, 

therefore, either U n or P„ must be specified along the outflow boundary, 

" oi a c 
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SPECIAL FORMS OF THE POISSON EQUATION FOR CONSTRAINT 
PRESSURE ALONG LATERAL BOUNDARIES 


The central difference approximation of the governing equation for the constraint 
pressure, equation (41), takes on a symmetrical form in the problem region interior. 
Near real and open lateral boundaries, this equation loses its symmetry and assumes 
special forms which are studied in this appendix. 

Equation (41), the governing relation for P c , is 


D f " Dt ,119 \/I T x3 Lt+rfi \ ^3 = ±_ _s_/Va_ d ^c\ 

4t 2 2t /A 9x“ Js aa JO ' 01 ' s/A 9x a y g Q!Q! f 9 x°y 


(Bl) 


Consider the staggered time grid lattice with U^ directed along positive dx*, U 2 
directed along positive dx 2 , and x* = I Ax* and x 2 = J Ax 2 as shown in sketch (a). 


o 


X 

oxoxo *- 

x 2 , U 2 , J 

X 

0 

1 x 1 , Up I 

Sketch (a) 

The values of £, Up and U 2 are known at time t on the circles, and the values of 
£, Up and U 2 at time t + r are known on the crosses. From equation (Bl) the values 
of P c at time t are sought on the circles. 
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Interior Solutions 

When all the grid points are inside lateral boundaries, the finite -difference approx- 
imation for equation (Bl) ^for simplicity only one vertical layer is considered and 



D f (I,J) - D^J) ^3 


4-r 2 


2 t 


,t+T, 


? l+ 7 (i+i,j) Ujd+i.j) - ? t+T d-i,j) Ujd-ijj) 


2 Ax ] 


Ax J 

2 T 


ot+T, 


u 2 (i,j+i) - ? t+T (i,J-i) u 2 (i,j-i) 

2 Ax 2 


D f (I+l,J) 

P C (I+2,J) - P C (I,J) 

D f (I-l,J) 

P C (I,J) - P C (I-2,J) 

2 Ax 1 

2 Ax 1 

2 Ax* 

2 Ax^ 

D f (I,J+l) 

P c (I,J+2) - P c (I,J) 

D f (I,J-l) 

P C (I,J) - P c (I,J-2) 

2 Ax 2 

2 Ax 2 ^ 

2 Ax 2 

2 Ax 2 

' — > 


Closed Lateral Boundaries 

In order to describe the treatment of computations near closed lateral boundaries, 
the boundary parallel to dx 2 is used as an example. The computations at other bound- 
aries are treated similarly. Assume that the zero depth point associated with a closed 
lateral boundary occurs at 1-2 (see sketch (a)). Then the boundary is parallel to 
dx 2 at I - 1, and to prevent mass loss from the system, U^ +T (I-1,J) must vanish. 

For the single layer being considered and with y/A = = \J%22 = e( l uat i on (38) inte- 

grated over the depth gives 



The finite -difference notation is ^since Uj t+T (I-l,J) = 0^ 

^? t+T (I-l,J) Ujd-M) = D f (I-l,J)(P c (I,J) - P c (I-2,jg (B4) 
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When equation (B4) is substituted into equation (B2), the terms involving U 1 (I-1,J) and 
P C (I,J) - P C (I-2,J) cancel leaving 


D f (I,J) - d\i,J) 
4 t 2 


Ax 2 

? t+T d+i,j) Ujd+1,0) 

2t 

^ 2 Axl 


D f (I+l,J) 

P C (I+2,J) - P C (I,J) 

Ax 3 

? t+T (i,j+i) u 2 (i,j+i) - c t+T (i,J-i) u 2 (i,J-i) 

2 Ax* 

2 Ax 1 

2r 

2 Ax 2 


D f (I,J+l) 

P c (I,J+2) - P C (I,J) 

D f (I,J-l) 

P C (I,J) - P c (I,J-2) 

2 Ax 2 

2 Ax 2 

2 Ax 2 

2 Ax 2 


(B5) 


which does not involve U-^(I-1,J) or P C (I-2,J). As an alternative to using equation (B5), 
one can use equation (B2) for the case under consideration by applying the fictitious 
boundary conditions 


U^I-M) = 0 ( 

P C (I-2,J) = P C (I,J) 


(B6) 


Inflow Boundaries 


To describe boundary computations for inflow boundaries, a boundary passing 
through the grid point (I- 1 , J) parallel to dx 2 is used as an example. The development 
of an equation for the computation of P c follows the development of equation (41) in the 
text. Consider a finite-difference approximation for equation (40) ^for one layer and 

\/A = Jg^ = Jgg 2 = l) which contains both boundary values and interior values of U a 


D f (I,J) - D^J) 


2r 


= -Ax' 


? t+T (I+l,J) U^+^I+^J) - ? t+T (I-l,J) U 1 t+T (I-l,J) 

2 Ax* 


- Ax' 


? t+T (I,J+l) U 2 t+T (I,J+l) - ? t+T (I,J-l) U 2 t+T (I,J-l) 

2 Ax 2 


(B7) 
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In the case under consideration, Uj^ +T (I-l,J) is a boundary value which, according to 
appendix A, must be prescribed. The other three velocities involved in equation (B7) are 

U 1 t+T (I+l,J), U2 * +t (I,J+ 1), and U2*" +T (I,J-1). Each of these velocities must be computed 
from either equation (38) or else equation (39). Recalling again the simplification used 

^one layer and \JK = = \J§22 = 1 )’ one has for u j t+T (I+l,J) from equation (38) 


Ax° 

2t 


? t+T d+i,j) u^i+i,^ = ^ ? t+T d+i,j) Ujd+ijj) - Df 2 ( -^ T~[ p c (I+2 ’ J) ■ p c (I ’ J 3 

(B8) 


for U2 ^ +t (I,J+ 1) from equation (39) 


^ C t+T (i,j+i) u 2 t+T (i, j+i) = ^ ? t+T d,J+i) u 2 d,J+i) - - p c (i,J)) 

(B9) 


and for U 2 * +T (I,J- 1 ) 
AX 3 lot+T/r t i\ tt t+r 


— ? t+T (I,J-l) U 2 t+T (I,J-l) = ^ ? t+T (I,J-l) U 2 (I,J-1) - Df(I,J l (p (l,j) _ P (I,J-2j) 
r z 2 t ^2 Ax' 2 


(BIO) 


Next, the governing equation for P c for the case under consideration is obtained by 
substituting equations (B8), (B9), and (BIO) into equation (B7) to get 


D f (I,J) - D t (I,J) ^3 

4^2 + ~27 


? t+T (I+l,J) Ujd+^J) - ? t+T (I-l,J) U^+^I-^J) 


2 Ax 1 


Ax c 


2t 


? t+T (I,J+l) U 2 (I,J+1) - e L+T (i,j-i) u 2 (I ’ J - 1) 


t+T/ 


2 Ax' 1 


D f (I+l,J) 

P C (I+2,J) - P C (I,J) 

D f (I, J+l) 

P c (I,J+2) - P C (I,J) 

2 Ax 1 

2 Ax 1 

2 Ax 2 

2 Ax 2 

D f (I,J-D 

P C (I,J) - P c (I,J-2) 



2 Ax 2 

2 Ax 2 




(Bll) 
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Equation (Bll) is the proper form of equation (41) to use adjacent to inflow boundaries. 
Here again one can make equation (B2) equal to equation (Bll) by specifying 

U^I-M) = U 1 t+T (I-l,J) (B12) 

along with 

P C (I-2,J) = P C (I,J) (B13) 

as a fictitious boundary condition on P c . Thereby the symmetry of equation (B2) is 
preserved. 

Outflow Boundaries 

The constraint pressure is usually specified along outflow boundaries. Velocity 
boundary conditions are then specified by upstream differencing. 
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Figure 2.- Current vector field from model in quasi-equilibrium. 10 days from initiation 
One (two) bars across arrow denote two (four) times scale. 
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